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ABSTRACT 


Several  methods  for  finding  eigenvalues  of  sym¬ 
metric  five-diagonal  matrices  are  compared  experi¬ 
mentally.  The  results  indicate  that  if  relatively  few 
eigenvalues  are  desired  a  modified  Sturm  sequence  and 
interpolation  scheme  is  fastest. 
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1.  INTRODUCTION 


We  are  concerned  here  with  finding  the  lower  (or 
higher)  eigenvalues  of  symmetric,  five-diagonai  matrices. 
Several  numerical  techniques  are  experimentally  com¬ 
pared,  including  one  that  uses  a  Sturm  sequence-type  ap¬ 
proach  directly  on  the  five-diagonal  matrix.  The  results 
indicate  that  if  relatively  few  eigenvalues  are  desired  a 
modified  Sturm  sequence  and  interpolation  scheme  is 
fastest. 
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2.  LR  TRANSFORMATION  OF  RUTISHAUSER  (LR) 


This  method  uses  Cholesky  decomposition  with 
appropriate  origin  shift  (Refs.  1  and  2,  and  Ref.  3,  pp. 
544-556).  Each  iteration  requires  6n  additions.  On  multi¬ 
plications,  2n  divisions,  and  n  square  roots.  The  method 
also  requires  a  judicious  choice  of  origin  shift  to  be  effec¬ 
tive.  A  poor  choice  could  invalidate  the  technique.  The 
implementation  used  here  was  that  of  Ref.  2,  which  pro¬ 
duces  a  cubically  convergent  technique.  Since  only  a  few 
eigenvalues  are  needed,  we  have  O(n)  operations.  Although 
the  Cholesky  decomposition  is  used,  the  origin  shift  allows 
one  to  find  eigenvalues  of  matrices  that  are  not  necessarily 
positive  definite. 
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3.  STURM  SEQUENCE  AND  BISECTION  (SS) 


By  the  separation  theorem,  it  is  known  that  the 
eigenvalues  of  the  (k-l)th  leading  principal  minor,  say 
A^_i,  of  a  symmetric  matrix  separate  the  eigenvalues  of 
the  kth  leading  principal  minor,  A^,  for  each  k  £  n  (see, 
e,  g.  ,  Ref.  3,  p.  103).  Thus,  if  we  evaluate  the  leading 
principal  minors  of  (A -XI),  the  number  of  variations  in 
sign  is  the  number  of  eigenvalues  of  A  greater  than  X. 


Let 


In  Ref.  4,  Kuttler  presented  the  following  recurrence  re¬ 
lation,  wh'  ~e  mi  is  the  ith  principal  minor: 


m 


k-1 


k  s  0 


m  =1 
o 


qk-2  =  bk-l  mk-3_Ck-l  \-3 


mk  =  akmk-rbkmk-2_Ck(ak-1  mk-3-Ck-l  mk-4>  +  2bkCkqk-2 


k  =  1,2, 


n  . 
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A  count  of  the  number  ot  operations  shows  that  we  have 
roughly  lOn  multiplications  and  5n  additions,  or  l5n  opera¬ 
tions  per  iteration.  It  is  necessary,  however,  to  con¬ 
stantly  scale  the  process  to  prevent  underflow  or  overflow 
(see  also  comments  at  end  of  Section  6). 

In  Ref.  5  Sweet  presents  another  set  of  recurrence 
relations  to  evaluate  successive  minor.  However,  sim¬ 
plifying  his  equations  for  a  symmetric  matrix,  it  turns  out 
that  one  needs  roughly  14n  multiplications,  8n  additions, 
and  n  divisions,  or  23n  operations  per  iteration,  which  is 
not  competitive  with  the  above. 
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4.  STURM  SEQUENCE  AND  INTERPOLATION  (SSI) 


In  Ref.  6,  a  scheme  is  described  whereby  the  num¬ 
ber  of  bisections  (and  functional  evaluations)  may  be  re¬ 
duced.  Once  bisection  has  isolated  a  single  eigenvalue  in 
some  interval,  an  interpolation  procedure  is  initiated. 

The  procedure  is  attributed  to  Wyngaarden,  Zonneveld, 
Dykstra,  and  Dekker,  and  its  implementation  is  described 
in  Ref.  6.  Since  the  scheme  requires  a  single  eigenvalue 
to  be  isolated,  one  cannot  find  multiple  eigenvalues  by  this 
method. 
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5.  BAND  DECOMPOSITION  AND  INTERPOLATION  (BD) 


In  Ref.  7,  a  scheme  is  given  whereby  a  band  matrix 
is  decomposed  into  its  LU  product  with  pivoting.  If  the 
matrix  is  symmetric  and  if  (A -XI)  is  decomposed,  one  can 
determine  the  variations  in  the  signs  of  the  principal 
minors  and  thus  use  the  approach  of  Section  4.  Since  the 
decomposition  requires  7n  additions,  7n  multiplications, 
and  2n  divisions  or  1  6  n  operations  for  a  five-diagonal 
matrix ,  the  method  appears  to  be  competitive  with  that  of 
Section  3.  We  have  not,  however,  taken  into  account  the 
bookkeeping  necessary  for  pivot  determination,  interchang¬ 
ing,  etc. 
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6.  MODIFIED  STURM  SEQUENCE  (MSS) 


In  an  effort  to  avoid  the  constant  scaling  of  the 
Sturm  sequence  relation,  we  consider  the  approach  used 
in  Ref.  8.  Let 


Mk 


m. 


m 


k-1 


Q, 


k-2 


\-2 


m 


k-1 


Then  Eq.  (1)  becomes 


M1  =ai 


M2  =  a2  "  M, 


Q,  = 


1  MxM2 


.  2  2 

M  b3+C3a2 

M3_a3  M2  MlM2 


+  2b3  C3  Q1 


Q 


k-1 


Vi 


k"2  Mk-1  Mk-2 


k-1  M 


k-1 


M.  =  a.  -  — 
k  k  M 


k-1 


Mk-1  Mk-2  \ 


+  2bkck«k. 
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However,  it  was  pointed  out  to  the  author  that  this  could  be 
greatly  simplified  as  follows: 


M1  =al 


P2=b2 


M2=a2-P2/Ml 


Pk  bk  “  Ck  Pk-l/Mk-2 


M.  =  a  -  c2/M,  _  -  P.2/M, 
k  k  k  k-2  k  k-1 


k  =  3, 


where  we  identify 

«k-i  ■  VVW  • 

This  is  equivalent  to  decomposing  the  matrix  into  its  LU 
product,  without  pivoting.  We  need  only  count  the  num¬ 
ber  of  Mj  <  0.  A  count  of  the  number  of  operations  in¬ 
volved  leads  to  2n  multiplications,  3n  additions,  and  3n 
divisions,  or  about  8n  operations  with  no  scaling.  Un¬ 
questionably,  this  is  the  fastest  method,  but  the  lack  of 
pivoting  introduces  a  question  of  accuracy. 

If  the  interpolation  method  of  Section  4  is  to  be 

n 

used  (SSI),  we  must  compute  v  M.  for  the  functional 

i=l  1 

values.  It  is  somewhat  simpler  to  scale  this  running 
produc  than  to  scale  the  recurrence  process  (Eq.  (1)). 

The  implementation  of  Eq.  (3)  was  done  in  a  man¬ 
ner  similar  to  that  described  in  Ref.  8.  A  test  was  made 
to  determine  if  vanished  and,  if  so,  to  replace  the 
zero  with  a  relatively  small  number.  This  never  occurred 


(3) 
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in  any  of  the  matrices  tested.  Also,  unlike  the  Sturm 
sequence  technique  for  tridiagonal  matrices,  it  is  not 
easy  to  see  the  conditions  under  which  two  successive 
(or  rrij  in  Eq.  (1))  vanish.  A  check  was  made  to  deter¬ 
mine  if  this  happened  and,  if  so,  to  perturb  X  and  try 
again.  However,  the  situation  never  arose. 
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7.  JACOBI  ROTATIONS  TO  TRIDIAGONAL  MATRIX  (RJ) 


If  all  eigenvalues  are  desired,  we  can,  of  course, 
use  any  of  the  above  methods.  However,  there  is  an  ap¬ 
proach  available  in  this  case  that  is  not  competitive  if 
only  a  few  eigenvalues  are  desired.  This  is  the  method 
of  Jacobi  rotations  to  reduce  the  bandwidth  of  a  matrix  but 
preserve  the  band  nature  of  the  matrix,  as  suggested  by 
Rutishauser  (Ref.  9)  (see  alto  Ref.  10),  followed  by  find¬ 
ing  the  eigenvalues  of  the  resulting  tridiagonal  matrix  by 
the  QR  method  (Ref.  11). 
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8.  RESULTS 


For  test  matrices  we  consider  the  following: 


/Vr 

2q 

r  \ 

2q 

p 

2q  r 

r 

2q 

p  2q  r 

A  = 

r 

i 

P  2q 

V 

r  2q  p-r/ 

/  n 

The  eigenvalues  of  this  matrix  are 

Xk  =  (p-2r)  -  i  (q2-(q-2r  cos  kfl)2) 

2  7T 

=  p  -  2r  -  4(q  cos  kfl  -  r  cos  kfl)  ,  8  = 

=  p  -  4q  cos  kfl  +  2r  cos  2  kfl  ,  k  =  l,  2, .  .  .  ,  n 

(see,  e.  g.  ,  Ref.  12).  In  particular  we  chose  four  matrices, 
as  follows: 


Matrix  1 : 

II 

a 

q  =  1 .  75  , 

r  =  0.  4 

Matrix  2: 

p  =  6  , 

q  =  1 .  75  , 

r  =  0.  5 

Matrix  3: 

P  =  11. 

q  =  10'15, 

r  =  5 

Matrix  4: 

p  =  10, 

q  =  IQ'15, 

r  =  5  . 

For  matrices  1  and  2,  the  eigenvalues  are  distinct 
but  of  the  smaller  ones,  many  are  near  1  for  matrix  1,  and 
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near  0  for  matrix  2.  For  matrices  3  and  4  each  root  is 
at  least  a  double  root  (for  even  n)  to  within  the  precision 
of  the  arithmetic  used,  with  distributions  similar  to 
matrices  1  and  2,  respectively. 

For  computing  only  a  few  eigenvalues  we  used  the 
methods  of  Sections  2,  5,  and  6.  For  computing  all 
eigenvalues,  we  compare  the  best  of  the  above  methods 
with  the  method  of  Section  7. 

Each  of  these  methods  was  applied  in  double  pre¬ 
cision  in  Fortran  IV  (H  compiler)  on  the  IBM  360/91  at 
The  Johns  Hopkins  University  Applied  Physics  Laboratory. 
Tables  1  and  2  contain  the  central  processing  unit  (C°U) 
time  in  seconds  for  each  run.  Because  of  the  time-snaring 
capabilities  of  the  machine,  it  is  possible  that  two  identi¬ 
cal  runs  could  show  slightly  different  times.  However, 
the  relative  magnitudes  of  the  numbers  are  consistent, 
and  it  is  these  in  which  we  are  interested. 

In  Table  1,  the  time  is  given  in  which  the  ten 
smallest  eigenvalues  were  found.  As  noted  above,  methods 
BD  and  MSSI  are  identical  except  that  the  Sturm  sequence 
process  is  determined  by  LU  decomposition  with  pivoting 
in  one  case  (BD),  and  without  pivoting  in  the  other  (MSSI). 
Operation  count  indicates  that  BD  should  take  about  twice 
as  long  as  MSSI,  but  the  actual  figure  was  between  six 
and  seven  times  as  long.  This  was  probably  because  of 
the  work  needed  to  determine  and  execute  the  pivoting. 

Also  examination  of  the  eigenvalues  of  the  matrices  up  to 
order  2000  showed  no  difference  in  the  accuracy  of  the 
two  methods. 

In  Table  2,  we  have  the  time  in  seconds  to  compute 
all  eigenvalues.  For  a  given  tridiagonal  matrix,  the  QR 
method  is  fastest.  We  have  included  the  QD  method  (Ref. 
13)  for  comparison.  It  would  seem  that  the  QR  method 
could  be  used  even  when  only  a  few  eigenvalues  arc  needed. 
Unfortunately,  "the  order  in  which  the  eigenvalues  are 
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Table  1 


r 

r . . 

Matrix 

rpr 

Time  in  Seconds  to  Compute  Ten 
Smallest  Eigenvalues 

i 

2 

i 

4 

Bf> 

6.  772 

6.  823 

9.  325 

3.  505 

N  -  500 

I.H 

2.  877 

2.  337 

1.  986 

1 . 463 

MSS 

2.  502 

2.  500 

1 . 374 

1 .  384 

i 

MSSI 

1.134 

1.195 

1 . 364 

1. 383 

i 

HI) 

1  3.  655 

1  3.  760 

1  7.  955 

18.  402 

|  n  -  i  ono 

I.H 

5.  887 

4.  629 

4.  196 

3.  01  3 

i 

MSS 

4.  851 

4.  91  3 

2.  610 

2.  760 

1 

MSSI 

2.  230 

2.  252 

2.  675 

2.  750 

i 

I.H 

12.  458 

5.  673 

8.  491 

5.  78.5 

x  H  2000 

MSS 

9.  223 

9.519 

5.  139 

5.  29  8 

MSSI 

4.412 

4.937 

5.  071 

|  N  =  5000 

MSSI 

10.  751 

10.  81  7 

11.  773 

■ 

1  2.  087 
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Table  2 


CPU  Time  in  Seconds  to  Compute  all  Eigenvalues 

_ _ _ 

Matrix  1 

Matrix  4 

N 

100 

200 

500 

100 

200 

500 

LR 

1.  369 

7.  278 

l.  311 

4.  935 

MSSI 

1. 835 

7.  277 

2.  818 

1  1.  108 

RJ  only 

0. 156 

0.  627 

3.  969 

0.  156 

0.  636 

4.  081 

RJ  -  QO 

0.  629 

2.  673 

1  7.  097 

1.  259 

5.  086 

RJ  +  QR 

0.445 

_ 

1.693 

.  .  .  .  ...  ^ 

10.  2.4 

0.  381 

1 .  523 

9.  348 
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found  is  to  some  extent  arbitrary"  (Ref.  11).  The  QD 
method,  however,  has  the  advantage  of  producing  the 
eigenvalues  in  strictly  ascending  order.  A  comparison 
of  Tables  1  and  2  shows  that  the  QD  method  might  be  com¬ 
petitive  if  more  than  a  few,  but  not  all,  eigenvalues  are 
desired. 


-  15  - 


TMt  JOHNS  HOPKINS  UNIVIASlTV 

APPLIED  PHYSICS  LABORATORY 

SlLVtN  IMUNO  MAKVLANO 


9.  CONCLUSIONS 


If  only  a  few  eigenvalues  of  a  large  five-diagonal 
symmetric  matrix  are  wanted,  the  modified  Sturm  sequence 
and  interpolation  technique  (MSSI)  appears  optimal.  If  all 
eigenvalues  are  desired,  Jacobi  rotations,  followed  by  the 
QR  method,  appear  difficult  to  beat.  Between  these  ex¬ 
tremes  on 'j  might  consider  the  LR  or  QD  methods. 
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